The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48_1/jurjen.fft

Internet Message Format | 1995-03-23 | 4KB

From en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!sdd.hp.com!spool.mu.edu!uunet!mcsun!hp4nl!charon!jurjen 18 Feb 91 10:28:30 GMT
Path: en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!sdd.hp.com!spool.mu.edu!uunet!mcsun!hp4nl!charon!jurjen
From: [e-mail] (Jurjen NE Bos)
Newsgroups: comp.sys.handhelds
Subject: HP28/48: Very FFT
Keywords: FFT
Message-ID: <[e-mail]>
Date: 18 Feb 91 10:28:30 GMT
Sender: [e-mail]
Organization: STORC, Veldhoven
Lines: 83
Originator: [e-mail]

A few friends asked me to write an FFT (fast fourier transform; converts a
vector to its "frequency diagram") program.  Looking through my archives,
I could not find any, although I vaguely remember another program was posted.
This program is written with speed in mind.  Most computations are done with
vectors instead of numbers.
At the end of the posting you see a HP48 directory.  HP28 users must know the
following:
The format is DIR <name> <value> <name> <value> ... END.
It is sent in translate code 3; that means that \.. codes are to be
translated:
\<< open program;        \>> close program;
\|v down arrow (make it a 'd')    \Gr greek rho (make it an 'R')
\GS greek Sigma (don't type in the routine DFT)
@ signifies comment.  Don't type in from here to end of line

A very nice feature of this program is that it works for ANY vector length.
If the vector length is odd, it is only slightly faster than the regular DFT
program in the end of the directory.  If the length is even, it is slightly
faster, and the more factor of two, the faster.  If the length is a power of
two, the regular FFT algorithm is applied.
All this is in one algorithm, without discrimination of all cases!

%%HP: T(3)F(,);
DIR
  CST { { "FFT"
    \<< 1 FFT
    \>> } { "-FFT"
    \<< -1 FFT
    \>> } }
  FFT    @ FFT routine.
    @ Input:    2: Vector 1: 1 for FFT, -1 for inverse
    @ Ouput:    1: transformed vector
    \<< OVER SIZE 1 GET
      IF OVER 0 <    @ Divide by length if inverse FFT
      THEN ROT OVER / ROT ROT
      END (0;6,28318530718) ROT * OVER / SWAP DUP
      WHILE DUP 2 MOD NOT    @ Divide out factors of 2
      REPEAT 2 /
      END SWAP OVER / \-> \Gr r p    @ p: twopower, r: odd rest
                    @ rho: root of unity
      \<< { r p } RDM
        IF r 1 \=/    @ Compute regular r-point FFT p times simultaneously
        THEN A\->L \Gr p * EXP 1 \-> m \Grp \Grn
          \<< 1 r
            START m 1 GET 1 2 r
              FOR k \Grn * m k GET OVER * ROT + SWAP
              NEXT DROP OBJ\-> DROP \Grp '\Grn' STO*
            NEXT
          \>> { r p } \->ARRY
        END TRN CONJ    @ TRN does transpose & conjugate, so conjugate back
        WHILE p 1 \=/    @ combine r-point FFTs to 2r-point FFTs
        REPEAT OBJ\-> DROP 'p' 2 STO/ { p r } \->ARRY TRN A\->L
    @ from here on, we work with conjugated numbers!
      \Gr NEG p * EXP 1 \-> m \Grp \Grn
          \<< { p r } \->ARRY TRN 1 r
            FOR k m k GET \Grn * OBJ\-> DROP \Grp '\Grn' STO*
            NEXT
          \>> { r p } \->ARRY + LASTARG - \-> m
          \<< OBJ\-> DROP m
          \>> OBJ\-> DROP 'r' 2 STO* { r p } \->ARRY TRN
      @ numbers are normal again
        END { r } RDM
      \>>
    \>>
  A\->L    @ convert matrix to list of rows.  This doesn't use sigma+, because
      @ this trick doesn't work with complex numbers :-(
    \<< OBJ\-> OBJ\-> DROP { } \-> r c u
      \<< 1 r
        START c \->ARRY 'u' STO+
        NEXT u
      \>>
    \>>
  DFT    @ The regular Fourier Transform.  Extra short version.
    \<< SWAP DUP SIZE 1 GET (0;6,28318530718) OVER / 4 PICK * \-> d x s q
      \<< 0 s 1 -
        FOR k '\GS(n=0;s-1;x(n+1)*EXP(q*k*n))' \->NUM
        NEXT s \->ARRY
        IF d -1 ==
        THEN s /
        END
      \>>
    \>>
END